function   modelOCP ( endTime )
% Get information of the OCP according to case$x.m
%
global sbspara NTswitch                 % by case$x.m
global nB nG iLD iGB_cut iG_cut idx_gx idx_gy  % by initial.m 
global omgs Pm EXd TJ TJsum             % by initial.m
global Pgen_cut Pld_cut                 % by initial.m
global Y2 YxV YxVfun Y3                 % by Ymatrix.m
global dataSim dlt_cut omg_cut V_cut    % read data or by simulate.m

global uGall uGalla I Ia                % global vector to reuse
global prob;                            % OCP problem
global uGez
% necessary parameters
dlt_max = sbspara(2);
I = zeros(2*nB, 1); uGall = zeros(nG, 1);
load(dataSim, 'dlt_cut', 'omg_cut', 'V_cut');
%% problem
prob.nG = nG; prob.nB = nB;
prob.ncutG = length(iGB_cut); prob.ncutL = length(iLD);

prob.time = [NTswitch(6) endTime];      % End time
prob.nstate = 2 * nG;                   % delta(1:nG), omega(1:nG)
prob.nctrl = 2 * nB + 1;                % Vx(1:nB), Vy(1:nB), COI(dlt)
prob.nparam = prob.ncutL + prob.ncutG;
prob.nocpvar = prob.nstate + prob.nctrl + prob.nparam;
% path: 0 = I - Y*V; 0 = coi - COI(delta); | delta - coi | <= dlt_max
prob.npath = 2 * nB + 1 + nG;
prob.ndae = prob.nstate + prob.npath;
% box bounds: {x,u}(T0) = {x,u}_T0; 0 <= uL, uG <= 1
prob.bnds.state_0 = repmat([dlt_cut; omg_cut], 1,2);
prob.bnds.ctrl_0 = repmat([V_cut; TJ' * dlt_cut/TJsum], 1,2);
%prob.bnds.ctrl_0 = [repmat(V_cut, 1,2); [-inf, inf]];
uGmax = prob.uGmax; uLmax = 1;                 % 0 <= uG <= uGmax; 0 <= uL <= 1
uGu = uGmax * ones(prob.ncutG, 1); uGu(uGez) = 0.001;
prob.bnds.params = [[zeros(prob.ncutG, 1), uGu]; ...
    repmat([0 uLmax], prob.ncutL, 1)];
% path bounds: -dlt_max <= delta - coi <= dlt_max
prob.bnds.paths = [repmat(0*[-1,1], 2*nB+1, 1)%zeros(2*nB + 1, 2); ...
    repmat(dlt_max*[-1, 1], nG,1)];
% objective & DAE
prob.objfun = @Obj;
prob.daefun = @Dae;

%% functions for objective and constrains
    function f = Obj (uG, uL)
        f = uG' * Pgen_cut + uL' * Pld_cut; % change to MW
    end
    function dae = Dae(dlt, omg, V, coi, uG, uL)
        Pe = EXd.*(V(idx_gx).*sin(dlt) - V(idx_gy).*cos(dlt));     
        if isa(V, 'double') == 1
            uGall(iG_cut) = uG; uEXd = (1 - uGall) .* EXd;
            I(idx_gx) = uEXd .* sin(dlt); I(idx_gy) = -uEXd .* cos(dlt);
            net = I - Y3 * V;
        else                            % !!! Ia for "auto-diff"
            feval(YxVfun, V, uG, uL);
            uGalla(iG_cut) = uG; uEXd = (1 - uGalla) .* EXd;
            Ia(idx_gx) = uEXd.*sin(dlt); Ia(idx_gy) = -uEXd.*cos(dlt);
            net = Ia - Y2 * V + YxV;    % for auto-differentiation
        end
        dae = [omgs*(omg - 1);       ...  % nG      ODE : 0 ~ 0
            (Pm - Pe)./TJ;           ...  % nG      ODE : 0 ~ 0
            net;                     ...  % 2 * nB  path: 0 ~ 0
            TJsum * coi - TJ' * dlt; ...  % 1       path: 0 ~ 0
            coi - dlt  ];               % nG      path: -dlt_max ~ dlt_max
    end
end